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Abstract 

Pulsc-pilcup affects most photon counting systems and occurs when photon detections occur faster than 
the detector's registration and recovery time. At high input rates, shaped pulses interfere and the source 
spectrum, as well as intensity information, get distorted. For instruments using bipolar pulse shaping 
there are two aspects to consider: 'peak' and 'tail' pileup effects, which raise and lower the measured energy, 
respectively. Peak effects have been extensively modeled in the past. Tail effects have garnered less attention 
due to the increased complexity: bipolar tails mean the tail pulse-height measurement depends on events in 
more than one time interval. We leverage previous work to derive an accurate, semi-analytical prediction 
for peak and tail pileup, up to high orders. We use the true pulse shape from the detectors of the Fermi 
Gamma-ray Burst Monitor. The measured spectrum is calculated by writing exposure time as a state-space 
expansion of overlapping pileup states and is valid up to very high rates. This expansion models losses due 
to fixed and extendable deadtime by averaging overlap configurations. Additionally, the model correctly 
predicts energy-dependent losses due to tail subtraction (sub-threshold) effects. We discuss pileup losses in 
terms of the true rate of photon detections versus the recorded count rate. 

Keywords: pulse pileup, deadtime, paralyzable, GBM, TGFs, model 



1. Introduction 

Pulse pileup affects x- and gamma-ray counting 
systems in the presence of high-intensity sources. 
At high rates the detected events become increas- 
ingly different from those registered by the instru- 
ment Q Detected events produce shaped electronic 
pulses that require some minimum time interval to 
measure the pulse height, and additional time to 
recover. The total time elapsed is the 'pulse win- 
dow'. Additional detections in this interval will 
modify the expected pulse shape and distort the in- 
ferred pulse-height information and detection rate. 
In this paper we elaborate on previous methods 
to develop a correction for pileup using the true 
bipolar pulse shape and with fewer assumptions [T] . 



1 Terms such as detected, input, incident will refer to the 
photons initially interacting in the detector volume. Terms 
such as registered, recorded, measured, observed will refer to 
the subset of photons and pileup events which are finally 
counted by the instrument 



We apply this method to the space-borne gamma- 
ray instrument Fermi Gamma-ray Burst Monitor 
(GBM). GBM consists of a set of inorganic scintil- 
lation detectors and is one of two instruments on the 
Fermi gamma-ray space telescope. An overview of 
the GBM detector hardware, science mission, and 
the context for pulse-pileup studies is described in 
section 

In bipolar-pulsed instruments pileup has two ef- 
fects on the measured pulse spectrum, depending on 
the arrangement of pileup events. Peak effects arise 
from the addition of positive pulse sections, caus- 
ing a single apparent high-energy count. Tail ef- 
fects occur when events are detected in the negative 
tail of the previous event, causing a second count 
with reduced apparent energy. For pulse shapes 
with large negative swings, modeling the tail pileup 
is especially important since losses occur when the 
summed peak is below threshold. Such losses de- 
pend on the input spectrum, a dependence which 
gets stronger as the bipolar peak-to-peak ratio ap- 
proaches unity. 
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With some assumptions the distorted spectrum 
can be predicted in terms of the true input rate and 
spectrum. Cano-Ott [2] and Danon [3] model peak 
effects in the context of unipolar pulses. Taguchi et 
al. model both the peak and tail effects in an X-ray 
counter with a bipolar pulse [T] . The approach gen- 
erally taken is to model the 'pileup response of the 
detector, which is a conditional probability (a.k.a. 
'likelihood function') of recording a certain energy 
£ resulting pileup of two events with energies E 
and E\: 

¥i{e\E ,E 1 ) 

Calculating such functions requires basic informa- 
tion about the hardware response of the detector to 
photon energy deposition, and a statistical model 
for the time-separation between detections. 

Assuming linearity in photo-current production 
and pulse shaping, the detector's pileup behavior 
is prescribed using a superposition of individual 
pulses. Then the likelihood functions are calculated 
assuming Poisson distributed events. The pileup 
spectrum is written using 'total probability', an in- 
tegral or sum of the pileup response into an mea- 
surement channel times the probability of input en- 
ergies (i.e., the input spectrum). Sections [6] and [7] 
derive the likelihood functions we use for modeling 
GBM pileup. 

A common simplification is to forego the true 
pulse shape when calculating the pileup-energy like- 
lihood, and instead use a mathematical approxima- 
tion. This has the advantage that likelihoods can 
have closed form expressions. For example, Taguchi 
et al. use a triangular or 'delta' approximation, 
so that the superimposed peak-time and amplitude 
are known and the dependence on event separa- 
tion is algebraically invertible [T]. In section 6.2.3 



we demonstrate why invertibility is necessary for a 
closed form likelihood expression. The disadvan- 
tage in approximating is that the true pulse shape 
gives significantly more accurate predictions than 
the approximate shapes. [2] demonstrate this by 
evaluating their model with an accurate gaussian 
pulse and then with several approximate shapes. 
They compare results with monte carlo simulation 
and show that only the true pulse shape is accurate 
over the entire energy range. However, they include 
only first order peak effects and reach a maximum 
pileup fraction well below what is expected in our 
applications. 

The technique of [1] allows for modeling high or- 
der pileup effects using an iterative approximation. 



Tail modeling uses convolution method rather than 
the pileup-likelihood technique used for the peak 
effect. This admits a number of simplifying as- 
sumptions about the timing and energy distribu- 
tion of tail events and avoids complications in cal- 
culating additional likelihood functions. However, 
when applied to the true pulse shape and specific 
input spectral shapes in our application, we found 
such assumptions did not yield sufficiently accurate 
results. This is perhaps due to a larger negative 
amplitude in our instrument, making the distorted 
spectrum more affected by tail pileup. 

The peak pileup method described in this paper 
is largely an application of the method of [T] to the 
GBM pulse shape. We present a different method 
for the tail modeling which extends the likelihood 
method to tail pileup, allowing us to model the 
spectral and timing randomness of tail events. This 
method also accounts for events occurring in a fixed 
deadtime that is separate from the peak pileup in- 
terval. 

A novelty of this approach is that the total spec- 
trum is obtained by partitioning time into overlap- 
ping pulse windows. Pulse windows are written as 
states of a Poisson process, each one having an asso- 
ciated peak and tail spectrum. In section [4] a state 
is defined by the number of events in the window, 
with pileup states having a non-zero number (the 
pileup order). In sections [6] and [7j the recorded 
energy distributions are derived based on the num- 
ber of events per state and how they are arranged 
into three sub-intervals of the pulse window. The 
total time in which there is a non-zero analog sig- 
nal is then the union of all pulse windows, and the 
total spectrum is the corresponding union of peak 
and tail spectra per state. This is expressed as su- 
perposition of pulse states minus overlapping terms 
(section [9]). We compare the model to monte carlo 
simulations in section ri0.3l 

As a final motivation for using the true bipolar 
pulse shape, we note there is only a small increase in 
computational overhead. Despite the apparent con- 
venience of closed form likelihood functions, the ac- 
tual difference in utility between them and a numer- 
ical function is small. Most spectra are sufficiently 
complicated that the total probability of likelihood 
times spectrum over energy is not analytically inte- 
grable, thus numerical evaluation is required. This 
is especially the case for non-ideal detectors, which 
have a complicated instrument response R(E,E y ), 
and detected counts spectra are already expressed 
numerically. In this case, other than some stor- 
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age and memory access overhead, computation of 
the pileup spectra requires about the same number 
of operations in either approach. As will be shown, 
for the constant Poisson process the likelihood func- 
tions are independent of source intensity and thus 
can be calculated a priori, kept in disk memory, and 
read by a separate program as required. In section 
|10| we describe a simple numerical calculation of the 
likelihood functions. 

2. Fermi-GBM 

GBM is a gamma-ray counting instrument 
aboard the Fermi Gamma-ray Space Telescope. It 
was launched in June 2008 and has been in near- 
continuos operation since the start of normal opera- 
tions one month later. GBM consists of 12 sodium- 
iodide (Nal) and two bismuth-germanate (BGO) 
detectors producing time and energy resolved data 
sets. Nal detectors have an effective energy range 
from approximately 8 keV to 1 MeV, and BGOs 
from 200 keV to 40 MeV. Pulse heights are digi- 
tized into 128 psuedo-logarithmic spectral channels 
per detector. More information on GBM hardware, 
detectors, and electronics is available in [4|. In- 
strument calibration (the relationship between peak 
voltage and input energy) studies are described in 

0-0 

Each detector has its own shaping and digitiza- 
tion firmware for performing pulse height analysis 
of detected events. Pulse-pileup, when it occurs, is 
on a detector-by-detector basis rather than an in- 
tegrated signal. GBM is not optimized for observa- 
tion of intense sources. The pulse shape has a finite 
width that requires 0.8/xs to register a single event 
and an additional 3.5 — for baseline recovery 
(section [3|. Generally speaking, pulse pileup oc- 
curs when the separation between detected events 
is smaller than the registration + recovery time. 
The negative peak amplitude is about 70% of the 
positive peak amplitude. To help mitigate pileup ef- 
fects, a fixed deadtime t is applied once a peak has 
been found, preventing pulse measurement during 
recovery. Pulse pileup however extends the time for 
baseline recovery and r becomes insufficient, adding 
tail distortions to the measured spectrum. 



2 For simplicity the figures in this paper are calculated 
using BGO channel energies that are approximated to be 
logarithmic: E m = (150keV)(1.04557) m gives the start of 
the m th energy channel in keV. The analogous Nal energy is 
given by E m = (4keV)* (1.04408)™. Actual channel energies 
are not precisely logarithmic. 



GBM is primarily an astrophysics experiment de- 
signed to detect transient cosmological sources such 
as gamma-ray bursts (GRBs). Detectors are uncol- 
limated and typical source + background rates are 
1 — 2 x 10 3 cps per detector. In the majority of 
cases a simple non-paralyzeble deadtime correction 
is sufficient, but pileup is potentially problematic 
for several types of sources observed. Depending on 
the application, pileup effects would become non- 
neglible at detection rates (Ao) between (roughly) 
5 x 10 4 — 10 5 cps in one detector, corresponding 
to observed counting rates (\R ec ) between 30,000 - 
70,000 cps in that detector. 

In Nal detectors the observed count rates from 
soft-gamma repeaters (SGRs) can peak at Xn ec ~ 
10 5 cps [5J. Solar flares are also routinely observed 
at these levels. Terrestrial gamma-ray flashes 
(TGFs) typically peak around Xr £C ~ 10 5 cps in 
BGO detectors [7]. At these levels the true de- 
tection rate can be 2 — 3 times higher, implying 
a substantial fraction (~65% - 75%) of the input 
photons pileup. Modeling the spectrum and ener- 
getics of such high-intensity transients thus requires 
understanding pileup distortion in GBM, and our 
goal is to have a fast, flexible numerical prediction 
which might be employed in spectral analysis. 

The model we present gives accurate results up 
to very high values of the input rate Ao- Demon- 
strations of spectral distortion versus rate are given 
in section [lOl The relation between \n ec and Ao is 
considered in more detail in section 1111 

3. Pulse shape 

GBM has well-studied electronic and detector 
properties that allow us to model its behavior with 
mathematical expressions. Shaping circuits make 
pulses whose first zero-crossing is fixed relative to 
the pulse start, regardless of the total photocurrent. 
We assume the zeroth order shape can be written 
in a separable form where the shape /(f) is mod- 
ulated by a scalar that depends on the amount of 
input energy. Thus a given pulse can be written 
V(t) — vof(t), where vq is an energy dependent 
coefficient, and /(f) is the basic pulse shape. Us- 
ing data collected in pre-launch testing, we fit the 
function below to sampled analog signal forms, for 
a range of input energies 

/(f) = K( Cl * t a - c 2 * f^e" 7 ' (1) 

where c\, C2, a, /?, 7 are fitted constants that param- 
eterize the pulse shape. K \s& normalizing constant 
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such that f(t) is unity at its peak (section 6.1 ). 

The best-fit shape parameters are c\ = 26.0, C2 = 
31.0, a = 1.27,/? = 3.5,7 = 2.6, if w 0.41023, for 
the time t is in microseconds. The figures and spe- 
cific results in this paper refer to equation ([I]) with 
these parameters. However, as an argument for the 
generality of the derived technique, we have tried a 
wide range of parameter values and find the model, 
when compared to monte carlo simulations using 
the pulse same shape, to be equally accurate in all 
cases. 
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Figure 1: A single unit pulse showing the zeroth-order peak 
time tp, the applied deadtime r, and the partitions A, B, 
and C. The pulse "window" is At = ta + TB + tq = 4.5 fis. 



3.1. Partition of the pulse window 

The pulse window initiated by detection of an 
initial event will be partitioned into three regions, 
each corresponding to different distortion effects. 
From the onset of the zeroth pulse, which can be 
arbitrarily set to to =0, pileup will occur only if 
there is another detected photon in one of these 
sub-intervals. We label them 'A', 'B', and 'C, hav- 
ing widths ta,tb, and tq, respectively (figure[l]). A 
and C will also be called the peak and tail regions. 
In our model, the particular distortion observed will 
be parameterized in terms of the number of events 
in each sub-interval. Choice of widths depends on 
the specifics of the pulse shape, the algorithm which 
measures pulse height, and the deadtime implemen- 
tation. However, we can give the following generic 



definitions for each region, based on the case of first- 
order pileup. 

Suppose we detect two photons within a pulse- 
width of each other. The initial pulse we call the 
'zeroth' event, and let it begin at time to = 0. Then 
the next event, beginning at t\ and having an as- 
sociated voltage V\, can occur in one of the three 
intervals A, B, or C. If t\ is in A, peaks will add and 
a 'summed' pulse height appears which is a function 
of Vq, Vi, and t\ — to- If ti is in B, vq is measured cor- 
rectly, and the second the peak is shifted down due 
to the negative baseline. In our case this coincides 
with deadtime, so there is no second count. Finally 
if t\ is in C, there is a second measurement, but 
the peak is shifted down from its nominal value V\ . 
Thus the widths are defined as follows, 'ta' is the 
minimum time required to resolve the pulse height 
of one event. For our model to apply it should be 
(approximately) < the positive signal width, 'tb' 
is a dead interval: if the second event occurs in this 
region it is not measured, but will influence the re- 
gion C. 're' is a live region, but lasts until the zeroth 
pulse is baseline-recovered. 

In the specific context of GBM, the width ta of 
the first region is the time-to-max tjj plus a buffer 
following the peak: 



T A — t p + TBuff 



0.78 fis 



This buffer is part of the instrument's digital peak 
finding algorithm, which requires a decreasing volt- 
age for a sequence of digital samples before regis- 
tering a pulse height j^j If additional events occur 
before the buffer expires (i.e., anywhere in region 
A), they will cause peak-pileup. [S] 

Events in part B are lost. Applied deadtime r 
contains the B interval, but is smaller than r be- 
cause a pulse occurring on the periphery of r might 
carry into the live time and be measured. Likewise, 
a pulse occurring in the post-peak buffer would be 
measured. Thus tb is the deadtime r less this buffer 
on the left, and the rise time on the right: 

tb = t - (t° + T Bu ff) = t - t a = 1.82 us 

where t = 2.6 fis is the GBM deadtime. 



3 This scheme is designed to eliminate measurements due 
to electrical noise. GBM samples analog pulses with a period 
of 0.104 /is. The buffer is programmed to be four samples 
long, so TBuff = 0.416/US. The deadtime results from waiting 
an additional 21 samples after a peak is found and buffered, 
so r = (4 + 21) * 0.104 = 2.6. 
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Region C is cutoff at an arbitrary point where 
the baseline (due to the initial pulse) can be con- 
sidered to be recovered. But the instrument is also 
live here, and peak pileup from additional events is 
possible. Thus we also require it to be short enough 
that only a single peak can be recorded. For GBM 
modeling we use the value tq — 1.9 /is, to give a 
total pulse window of width 



t~a + i~b +t c = At = 4.5/xs 



(2) 



4. Pulse-pileup 



Having partitioned the pulse shape into three ad- 
jacent regions, we can now describe pulse-pileup in 
terms of the number of events in each one. Assum- 
ing a zeroth event to begin the window, one or more 
added events in region A constitutes peak pileup. 
One or more events occurring in region C cause tail 
pileup, or the tail effect. Additionally, pulses in B 
contribute to baseline reductions for the C events, 
and this is a more extreme category of tail effects. 
Together these produce artificial raising and lower- 
ing distortions in the measured spectrum, which is 
our goal to model. 

Assuming a zeroth event with the A-B-C parti- 
tions, the state or configuration of pulse-pileup can 
be represented using three non-negative integers 
a, 6, c which give the number of additional events 
in each interval (excluding event 0). We define the 
symbol (a, b, c) to describe the pileup state in the 
pulse window Ar. The total order of pileup is the 
number of additional events, such that order zero 
is a single count without pileup, first order is one 
added event in A, B, or C, etc. Thus the generic 
configuration (a, b, c) represents the case of an inci- 
dent event (the zeroth event) followed by (a + b + c) 
events within the window. In later sections we ad- 
dress the probability and independence of pileup 
states, but for now we need only the definition: 



Event (a, b, c) 

= a pulses in A, b in B, and c in C 

In general this represents the event of (a + b + c) th 
order pileup, where (a + b + c) > 0, and (0, 0, 0) 
is the event of a single recorded count, accurately 
measured. In this scheme the maximum number 
of measured counts is two (one from the peak, one 
from the tail). 



Figure [2] demonstrates the simplest case: first- 
order pileup. Figure [3] gives several examples 
of higher order pileup, with the pulse partitions 
shown. The possibility of a pulse in B occurs at 
second order and above. In the following sections 
we calculate pulse height distributions correspond- 
ing to each state (a, b, c), and an input energy spec- 
trum. The peak spectrum will depend only on a, 
but the tail spectrum depends on all three orders. 

One consequence of pulse pileup is that the dead- 
time is extended beyond its the nominal value r. 
For example, pileup in region A results in a summed 
peak at a time > the single-pulse peak time 
And since the deadtime is imposed from this peak 
it effectively migrates forward, contributing to par- 
alyzable deadtime 9J. Strictly speaking the pulse 
partitions should have a corresponding shift which 
depends on this recorded peak time. Instead we 
use a semi-paralyzable model where the summed 



peak time (section 10.1 ) can wander into region B, 
but the partitions and deadtime remain fixed. The 
result is shown to be accurate when compared to 
monte carlo simulations. 

5. Interval Distribution 

Within fixed intervals the interval distribution 
(i.e., the probability density (PDF) of separation 
time between events) can be calculated under the 
assumption of a fixed-rate Poisson process. This 
derivation is well-known, appearing in [lj among 
other places, but because of its centrality in the 
calculation of pileup likelihood we also present it. 
For the Poisson process of rate A, the PDF of the 
separation s between any event and the next is 



f s (s) = Ae 



-As 



(3) 



which is defined on s e (0, oo). Define the separa- 
tion of the i th event from the (i — l)*' 1 as 

Si = ti - ti_! (4) 

where ti is the time of event i. Now suppose a finite 
interval ranging from to At, with the zeroth event 
outside the interval at to — 0. Then the probability 
of exactly one event in the interval is given by the 
(joint) probability that s± < At and si + S2 > At, 
and the second condition is equivalent to S2 > At — 

81. 

Pr(si < At and (s 2 > At - si) ) 

si oo 

. f . (5a) 



Ae- ASl d Sl 



Xe~ AS2 dso 



o 



At — si 




Figure 2: The three cases of 'first order' pileup, (100), (010), and (001), showing the measured peak for two events of equal 
energy, and the dead time r as imposed by GBM hardware, (a) shows the peak effect, and (c) the tail effect. Panel (b) depicts 
a nominal case where one count is accurately measured and the next is lost. Typically this is not regarded as 'pulse pileup' as 
there is no associated spectral distortion of the pulse height (only the count rate), which can be corrected by simpler means. 




Figure 3: Higher order pileup examples, with the A-B-C partitions shown, (a) second order peak pileup. (b) third order 
pileup, with peak and and tail effects, (c) a third order case of the deadtime+tail effect. Recorded pulse height in C depends 
on the tail from pulses in A and B. 
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Differentiating in s± gives the joint probability den- 
sity f s (si , f), whose general form is / s (si , n) for 
the random variables si and n (number of events 
in At) 

1) - ^"Pr - Ae- AAt (6) 
asi 

The conditional probability density, f s \ n {s), is given 
by dividing the joint density f s (si, 1) by the prob- 
ability that n = 1, a procedure equivalent to nor- 
malizing over the interval [0, At]: 



fs\n=l(s) 



fs(si,l) 



Ae 



-A At 



AAte 



— AAt 



1 

At 



(7) 
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Figure 4: 2 nd order pileup in the offset interval [t,t + At], 
showing events to,ti, t2>t3' If we assume finite pulses occur 
at each time, events t\,t% are piled- up, with possible tail 
effects from to. si = + t — to, and the random variables 
s' 1 ,S2 are identically distributed. 



which is a constant, and is independent of A. [TU] 

For exactly two events in the interval, the condi- 
tion is si + S2 < At and Si + s 2 + s 3 > At (e.g. 
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figure [4]). This can also be written s 2 < At 
and S3 > At — (si + S2). Then 



si 



Pr(si < At for 2 in At) 

si 

= I Xe- Xsi d Sl 







At-(si + s 2 ) 



A 2 e" AA *( Sl At-^) 



(8) 



Assuming a definite state of n = 2 in At gives 
the second order conditional density: 

/.(*!, 2) 2 



/s|2(Sl) 



/ At / s K,2)d S ; 



At 2 



(At - Sl ) (9) 



Because in this process event separations are inde- 
pendent, identically distributed random variables, 
si,S2 follow the same probability density / s | 2 . For 
the general case of n events, the joint density func- 
tion f s {s\,n) is calculated as in equation Q and 
the conditional density f s \ n (s) just like ([9|, which 
can be shown by induction to be, for all n > 1, 
n 



fs\n{ s l) 



At" 



(At - s a 



,n-l 



(10) 



5.1. Offset Intervals 

The above density was derived assuming to 



0, but equation (10 1 turns out to have the same 
form regardless of the position of the previous event 
to before the interval. Consider the second-order 
integral in equation (|8j), but now let the interval 
begin at some arbitrary time t > and cover t to 
t + At. Let the zeroth event occur at an arbitrary 
time to < t. Now the 2 nd order pileup conditions, 
for ti,t2 in the interval, are: 

t - t < si < (t + At) - t 

si + s 2 <(t + At) - t 

si + s 2 + S3 > (* + At) - t 

An example is shown in figure Q. Using the sym- 
bols T A = t-t and T B = (t+At)-t for shorthand, 
the second-order joint probability is 

Sl 

Pr(si,2) = J \e- Xsi d Sl 



Ta 



T B -si 



00 



I Ae 2 ds 2 j /ve -us 3 

T B -( Sl +s 2 ) 



Xe- Xs3 ds, 



si 



■2T B ( Sl -T A ) (11 



Differentiating and normalizing gives the condi- 
tional density: 

2(T B ~ Sl ) _ f[ t + M-t ]-8 1 ) 



(12) 



2(At- [si - (t-t )]) 
(At) 2 



Taking just the portion of si that's in the normal- 
ization interval, = Si — (t — to), the above is 



/ S | 2 (si) 



2(At- 



At 2 



(13) 



for Si < At, which is exactly the result when t 
is at the start of the interval, equation This 
generalizes to the result 



fs\n(s) 



At" 



(At - s) 



n-l 



(14) 



for n events in an arbitrary interval [t,t + At], 
and < s < At. Using this result we can write 
the probability densities assuming a definite pileup 
state (a, b, c): 



fs\a{s) = 


(r>*->" 


(15) 


f s \b{s) = 




(16) 


fs\c(s) = 


(rlr (T °- r ' 


(17) 



6. Peak pileup effect 

The peak modeling technique is an application 
of the one presented in pQ. The method is to first 
derive expressions of the form Vx(e\E' , E"), which 
give the likelihood of recording energy e when two 
events are detected with energies E',E", within a 
time ta of each other. For first-order the model 
is accurate without any approximation. At higher 
orders an iterative approximation is used since the 
dimensionality (i.e. number of random variables) 
becomes large. Then the total probability of e is 
given in terms of the pileup likelihood and the input 
spectrum. 
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separation of the two events, s±. The maxima of 



Figure 5: A case of 1 order peak pileup, showing the mod- 
eled peak at . 



6.1. Zeroth order peak-time 

Without pileup, a single pulse at time has a 
voltage signal that can be expressed as 



V(t) = v f(t) 



(18) 



where f(t) is the normalized pulse shape, such that 
at its maximum, f(t) = 1. Therefore the peak volt- 
age is V = vq. The peak-time is given by the first 
solution of 

df(t) 



dt 







(19) 



and we will denote it by the symbol t p . Thus, 
/'(*£) = 
f(t° P ) = 1 
V(t° p ) = v 

6.2. First order peak pileup 

6.2.1. First order peak-time 

First order peak-pileup is depicted in figure [2ja) 
and figure [5] The apparent pulse is a linear combi- 
nation of unit pulses f(t), 

V(t) = «o/(*-*o)+«i/(*-*i) (20) 

where to, t\ are the incident event times, and vq, v\ 
are their peak amplitudes. For the constant rate 
process, the total time offset of the pair is irrele- 
vant, so we can let t a = 0. Then t\ is equal to the 



equation ( 20 ) are given by 



dt 



V(t) 



t'o- 



df(t) , df(t- Sl ) 



dt 



■v%- 



dt 



(21) 



We define the symbol i~ as the position of the 
measured first-order peak, such that it satisfies the 
above equation. As discussed in section [4] the ob 
served peak is shifted forward, so tp 
depends on the random variables si,vq,V\: 



> t^, and it 



*p(si,i>o,i>i) 



(22) 



Having a model for ii(si, vq, v%) is necessary so that 
the sum ( 20 ) can be evaluated for the observed first- 
order energy. Depending on the pulse shape f(t), 



equation (21 1 might have a closed form solution for 
t^. In previous work simplifications are employed 
such that tp is given by simple geometrical argu- 
ments. But in the present case, substitution of the 



pulse shape, equation ([!]), into equation (21 1 results 
in an expression that cannot be readily solved for 
tp. As a result we use an empirical model for the 
peak-time that is presented in section [l0.1| 

Note that for a given set of values (si,vq,vx) 



there may be two solutions to equation (21) (two 
maxima in the peak interval, for example in fig- 
ure |5|. Selecting the correct peak requires knowl- 
edge of the peak-finding algorithm of the detector 
system in question. The buffering scheme of the 
GBM pulse-height analyzers generally causes the 
last maximum to be the measured one. However if 
vq is much larger than v±, and si is in the buffer, 
tp —> tp because the falling derivative has a much 
larger (negative) value than the rising derivative in 
this region. When the two add, the net change is 
still negative, so the decreasing sample criterion is 
satisfied. Such discrete logic is expressed as piece- 
wise behavior in the model for tp. 

6.2.2. First order peak- energy 

Let us assume for now that we have a sufficient 
model for tJ(si,«o,t>i). Unlike f{t° p ) = 1, /(**) £ 1 
in general. The recorded pulse height is 

V(t 1 p)=v f(t 1 p )+v 1 f(t 1 p-s 1 ) (23) 

and the recorded first-order pileup energy e\ is 

ei=Z[vof(t 1 p )+v 1 f(t 1 p -8 1 )] (24) 
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where is the channel or voltage-to-energy con- 
version, determined by standard methods of instru- 
ment calibration. For modeling input spectra, we 
assume the inverse £ _1 exists and the coefficients 
of input pulses can be calculated as Vi — £ _1 [£y 
where Ei is energy deposited by a detected gamma- 
ray. Ei is of course converted from incident photon 
energy by the various physical processes of the de- 
tector, but this is extraneous to the current prob- 
lem. It is sufficient to say that is the recorded 
energy of the i th count in the absence of pileup ef- 
fects. 

We can now give the general definition of the im- 
portant function si, which gives the recorded 1st- 
ordcr energy for a given instance of the random 
variables (si , Eq , E% ) , 



(25) 



£1(31,^0,^1) 



Note that if £[i>] is approximately of the form y 
mx, then this s implifies to E a f{tp) + Exf(t p — si). 



In section 10.1 we give the expressions for tp and £1 
using the true GBM pulse shape, which is used to 
calculate figure [6j 

6.2.3. Probability distribution of Ei 

What we need is a model giving the observed 
distribution of the recorded energy £1 in the event 
of peak pileup. This means we need the dis- 
tribution of £1 over all possible realizations of 
{(si, Eq, Ei)}. Let us examine the distribution £1 
due to s\ only, holding Eq,Ei fixed at some arbi- 
trary values. This gives a conditionalized form of 
equation (25), £i(si | E ,Ei). 

Given the interval PDF / s |„ = i(si) of equation 
(16 1, the PDF of the dependent variable £1 can be 
found in the standard way [10] - For si in the peak 
region, £i(si | Eq,E\) is monotonic, and we can 
write its probability density as 

1 1 \ \ fs\n=l( s l) 

fpJaki ) = d£l(s) (26) 

I ds I s=s 1 

where £i(si | Eq,Ei) is written £i(si) for brevity. 
It is understood that tp and £1 are functions of the 
three random variables. 

We can write the probability of pulse pileup into 
an energy interval £i(si) to £i(si) + de\ as 

d[^ { l k {ei{si) I E ,E 1 )] =/£L(£i)d£i (27) 



Then we can write the following for the probability 
that £1 is in the discrete channel e to e + Ae, in the 
event that it is due to Eq,Ei pileup: 



e+Ae 



/s|ra=l( s i) 



1 as 1 s—s' 



(28) 



ds' (29) 



This can be converted to an integral over si using 
the positive-definite Jacobian determinant (proba- 
bility must be positive) 



(30) 



d(e') 




dsi(s) 


d( Sl ) 




ds 



giving 



s(e+Ae) 

fs\lW 
I deiQ) I 
s(e) ' ds ' s i 

s(e+Ae) 



/ s |i0i)ds'a 



dei(s) 



ds 



(31) 



Examples of this likelihood function are shown in 
figure [8] Evaluation of the above integration limits 
evidently requires knowledge of an inverse function 
giving the separation in terms of the summed peak 
and two energies. Therefore, if a closed form ex- 



pression for equation (31) is desired, the recorded 
energy £i(si | Eq,Ei) must be an analytical expres- 
sion and have an inverse, Si = £|f 1 (£, -Et)i -^i)- 

The simplification of [1 using triangular pulses 
allows such an inversion. In our case, the true pulse 
shape is too complicated to invert, so instead we 
leave the expression in its implicit form and evalu- 
ate numerically. In section [T0| we describe a simple 
discrete algorithm for doing this implicit evaluation. 

Regardless, we can write the first-order peak 
pileup contribution to the spectrum as 



P(i)(si)= J J Pr£L(ei \E< h E 1 ) (32) 


x S(E )S(E 1 )dE dE 1 



9 



SiCsi) with E Q = 3500 keV, E x = 3500 keV 



, 4000 



■ Model for peak energy 

■ Model for tail energy 



Monte carlo measuremnt with phase uncertainty 
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Figure 6: Plot of lst-order recorded energy e.\{s\, Eo, Ei) for Eq, E\ at 3500 keV, if the second event occurs at s\. The curves 
represent the model for measured energy using the function tp, which is different in the peak than in the tail. Shaded areas 
between points give the range of recorded energy due to the phasing of digital samples. Discontinuity at the end of the peak 
occurs when the separation is large enough to allow pulse distinction by the instrument, and Eq is measured correctly. The 
divisions A, B, C refer to the pulse partition of section [3. 1 1 The curves in B-C refer to measurement of a second event (tail 
effect), which is addressed in section [7] 



S{E) is the PDF of the energy spectrum initially 
observed in the detector, and is normalized to unity. 
S{E) is typically modeled in terms of the externally 
incident photon spectrum and the detector response 
function R(E\E~), for example as described in [IT] : 

S(E) = J R(E\E J )S 1 (E 1 )dE 1 (33) 

where 5 7 is the PDF of the photon spectrum. Since 
real detector systems are sensitive over a finite en- 
ergy range, we note that the limits of integration in 
equation ( [32] ) can be finite values E m i n — > E max as 
long as S(E) is normalized over this interval. The 
PDF of the spectral contribution is 

p (1> (ei) = ^P (1> (ei) (34) 

6.3. Higher order peak terms 

For n th order pileup we require an expression giv- 
ing the recorded energy within the peak interval. 
The n + 1 pulses will be superimposed and recorded 
as a single pulse height e n . This energy would gen- 
erally be dependent on the energies and separations 



of the piled-up events: 

n 
i=0 

=> £ri = £n{si, ■ ■ ■ , S n | Eo, Ei,E 2 , ■ ■ ■ , E n ) 

for Si = U- ti-x 

The other expressions of the previous section would 
generalize in a similar way, leading to the n th order 
peak correction 

P{n) (En) 

= J dE J dE r -- J dE n 

x Pr£L( £ ™ I E ,...,E n )S(E )...S(E n ) 

(35) 

These expressions are complicated and unwieldy, so 
instead we use the iterative approximation of pQ. 
This scheme approximates higher order variations 
due to the added random variables by using results 
from the previous order. 

To calculate the n th order correction, the previ- 
ous order term p/ n _i) is used: 
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Figure 7: P r p ea fc(£|-£'Oi E\), for several combinations of Eq,E\. Assuming two peak pileup events, these recorded energy 
distributions result from varying their separation. The primary peak occurs at e = Eq + E\, indicating this is the most 
frequent case, but the width of each distribution demonstrates that Eq + E\ is not a good approximation for e in general. A 
second peak, appearing in panels (c), (d), (g), and (h), indicates that one of the two events has a higher probability of being 
distinguished. This occurs when one energy is much larger than the other and the separation is at least . 



P (n){£n) 







X P(n-l>(-En-l)<5(I? n )dEn-ld-En 

(36) 



and the kernel Pr^™^ is evaluated using the n th 
order interval statistics: 



s(e+Ae) 



Pr Sfe^ I E n -i,E n ) = J f s[n (s')ds f 



»00 



(38) 



with the PDF 



P(n)^n) = -7— P(n){Zn) 
Cl£ n 



(37) 



We make the approximation that e n ~ 
s n (s n | E n -i,E n ) and has the same form as 
£i(s n | E n -i,E n ). f s \ n {s) is the interval dis- 
tribution for n events in the peak, equation 



(16) 



We can unify the 1 order correction with (6.3) 



by defining p(fl) (E) = S(E). Then equation (6.3) 
reduces to d32l for n = 1. 
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7. Tail effect 



tl « t° + s±. Thus equation (25) becomes 



The "tail effect" is spectral distortion caused by 
the bipolar pulse tail, which reduces pulse heights 
occurring before baseline recovery. The model of 
Q] presents a clear and relatively accurate model 
for tail effects in the output spectrum. However, 
in that case there is not the intervening deadtime 
region (interval 'B'). For modeling purposes this 
means tail subtraction effects are taken to depend 
only on the peak state. In our case we must consider 
the events in the peak and deadtime, since pulses 
in either region contribute to a negative tail when 
the instrument again becomes live. The method 
of [T] uses further simplifications by reducing the 
number of random variables: input tail events are 
assumed to be uniformly distributed in time and 
have the same energy, which is taken to be the 
mean energy from the modeled input spectrum; i.e., 
S(E) -> S(E - (E)). 

We present an alternative tail technique which 
deals with the intervening deadtime interval (region 
'B' in figures [l] and [3]) , and models energy and tim- 
ing variations of tail pileup events. This method is 
similar to that used for peak pielup in that pileup 



(39) 



likelihood functions of the form of equation (31 ) are 
derived for measurement in the tail region (region 
'C'). Measurement in this region is sensitive to the 
previous two, so the likelihood scheme is more com- 
plex than in the peak case. The number of events 
in A and B affect the C measurement because their 
negative tails combine (figure [3jc)) . Additionally 
peak pileup in C (figure |3^b)) can occur and must 
be modeled. We account for the various possibilities 
by isolating an 'A+C effect, which assumes counts 
only in A and C, and a 'B+C effect, which assumes 
no zeroth event and counts only in B and C. Each 
effect is calculated for first order, then the iterative 
treatment is applied. The spectrum of the total 
pulse configuration, with 'A+B+C dependence, is 
calculated in terms of the two effects. 



7.1. Tail energy 



In region A we employed the function 
tp(s\, vq, vi) to give the time of the first-order 
peak. For tail pileup we make the simplifying 
assumption that the location of the summed peak, 
were it to be measured in B or C, is approxi- 
mately at the peak of the second input pulse, i.e., 



£i(si, E , Ei), for tail energies 



Figure [6] shows that this model is reasonably ac- 
curate. Note that this simplification is not neces- 
sary; if tp for the tail regions is known, equation 
( 25 1 can be used exactly. 



7.2. Recorded spectrum of tail events 

The recorded tail spectrum depends on the to- 
tal pulse the state (a, b, c) . We define the following 
symbols that will refers to the spectrum of mea- 
sured pulse heights for a given pulse state: 







probability of recording 






energy e in the tail 




( E ) = 


PDF of recorded energy, 







The affect of events only in A upon measure- 
ments in C is first expressed in the (a, 0, c) com- 
ponents. For c = there is no tail count and 
thus no tail measurement, so such terms are zero. 
For a + c > 1 an iterative scheme is used. In 
general, a succession of terms will be used to go 
from first order tail pileup to generic higher orders: 

Q(0,0,1) Q(0,0,c) ~~ * Q(a,0,c) ~^ Q(a,b,c)- 

7.2.1. First tail effect, Q( a ,o.c)( £ ) 

We begin by calculating the likelihood of record- 
ing energy e in region C in the event of tail pileup. 
Here we assume a zeroth event to begin the pulse, 
and additional events in A and C. For first order, 
a — and c = 1. Figure [9] depicts this configura- 
tion. 

Recall that the probability density of a single tail 
event is, f s \ c (f[) = ^, where s[ = Si - (t a + t b ) 
(equation ( 17) , section 5.1 ). The full separation s\ 



is used to calculate the measured energy e(s\) using 
equation ( [39] ). Assuming E in ta, Ei in tq, the 
probability of tail pileup into the discrete energy 
bin [e, e + As] is then: 

s(e+Ae) 

Pr£J. c (e I E ,E 1 ) = J f slc=1 (^)ds\ (40) 
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Figure 8: Recorded energy distributions for the two tail effects. In each plot the lighter histogram plots Pr^ , c (e\Eq, Ei), 
which assumes Eq is the peak energy, and the separation si ranges from ta + t~b to ta + tb + tq • The darker histogram plots 
PT 1 B+c (e\Eo, Ei), which assumes Eq is in the dead region (B), and si varies from tb to Tg + tq. 
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Figure 9: (0, 0, 1} pileup event, t p w tp + si for the tail pulse. 
s' x is distributed according to equation (TTTb. 



and the 1 st order spectral component is 

oo oo 

Q<0,0,1>( E ) = / / Pr A+ C (£ I ^0,^l) (41) 


x S(E )S(E 1 )dE dE 1 

For (0, 0, c) states and c > 1, we again do an itera- 
tive approximation using the previous order. How- 
ever, there are now multiple pulses in the region 
C, and we approximate their measured energy as a 
case of (c— l) th order peak pileup. But now the pri- 
mary and additional event have a negative baseline, 
and so we approximate them using <f(o o,c— l) a s the 
primary peak distribution^ and S(o,o,i) as the spec- 
trum of additional events. Then we use the kernel 



4 Much like the a th order peak iteration uses p a —i(e) in 
the primary position. 
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Pr 



<c-l) 



from equation ( 38 ) . This is an approxima- 



tion since Pr^g^j, models pileup in the interval ta 
instead of tq< However multiple events in tq tend 
to be clustered due to the distribution JW , and so 
the small difference between interval sizes ta and 
tc is not a major source of error. Then the (0, 0, c) 
tail spectrum is 



For c > 1, 




Q{o.o,c) (e) 

00 00 






(42) 







x 9(o,o,c-i>(-®')9(o,o,i; 


(E")dE'dE" 



For (a, 0, 1) states we approximate the effect of 
the multiple orders in A by using p^ in the inte- 
gral, and using the Pi^+c kernel, since the interval 
distribution in C does not depend on the number of 
events in A (section |5~T] ). The spectrum measured 
in the tail for these configurations is 



For a > 0, c = 1, 



where the c order kernel is approximated as 

s(e+Ae) 

Prf +C (e I E a ,E c )= j f^s^ds', (45) 



s(e) 



and the recorded energy is equation ( 39 ) with 
E a ,E c replacing E ,Ei and f s \ c is the tail inter- 



val density, equation (17) 
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Figure 10: Configuration used to calculate Pr^i. c (s\Eq, Ei), 



which assumes no zeroth event, ti. 



s' is distributed according to equation \17\ 



tp +si for the tail pulse. 



Q< a ,o,i>(e)= / / Pr$- C (e I E',E") 



(43) 







P{a) (E')S(E")dE'dE" 



For a = this reduces to equation (41) since 
P( }(E) = S(E). For zero events in the tail, there 
can be no recorded tail energy, so q/ a ,b,o}(E) = 
for all a, b. 

Finally, for (a, 0, c) states, we use the a th order 
peak spectrum with the ?(o,o,c) ( e ) just calculated: 

For a > 0, c> 1, 

00 00 

Q<.,o,c>(e) = J J Pr<g. c (e I 


xp (Q > (^')9(o,o,c-i> (E")dE'dE" 

(44) 



7.^.^. Second tail effect, Qt a ,b,c}( s ) 

The second effect models the case when events in 
B pileup and create a large negative peak in C. It 
is only present at orders > 2 since it requires events 
in B and C. We approximate the variations due to b 
random events in B using the peak spectrum terms 
calculated for order b — 1 (there is no zeroth event 
in B). Figure 10 depicts the lowest order B+C con- 



figuration, used to calculate Pr^+c varying s[ 



Equation ( 39 ) is used for the recorded energy, with 
the separation given as s — s' + tb, with s' dis- 



tributed according to equation (17). The probabil- 
ity of the pulse height lowering into [e, e + Ae], due 
just to an event in B, is 

Pr sic( £ I E b _i,E c ) = J f slc (s[)ds' (46) 

for < s < tc- 

Finally, the additional event in C is an event from 
the input spectrum, shifted down due to tails from 
A, but also shifted up due to pulses in C. In other 
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words, it is approximately distributed as 3( a ,o,c}( £ )- 
Thus we use these two terms with the kernel above 
to calculate the second component 

For b > 0, 

OO OO 

Q(a, b ,c)(s)= [ f P^ +C (S | E',E") 



xp {b -i ) (E')q {ai0te) (E")dE'dE" 

(47) 



This spectral component clearly depends on the 
total positive signal from both A and B. Thus we 
have an approximation for the highest order states 
using the succession of terms, Q(o,o,i) Q(o,o,c) 

Q(a,0,c) Q(a,b,c) ■ 

8. Probability of pileup events 

Since the time intervals A, B, and C do not over- 
lap, the number of events in each interval is inde- 
pendent of the other two. Furthermore, they are 
Poisson distributed random variables. We assume 
A is constant throughout the detection of the input 
spectrum S(E). The probability of each state is 
then 

Pr(state \ rate = A) 

= Pr(a G A and b £ B and c G C \ X) 

= (\T A n\T B )\\T C y c „ Mta+tb+tc) 

al *6! *c! 



= Pr«o,6,c)|A) 



(48) 



The rate A is separately modeled so the condition 
can be ignored and the expression for state proba- 
bility is 



Pr«o,M)= (WW e -^ 



a! *&! *c! 
for t a + t b + t c = At 



(49) 



This gives the probability of having non- 
overlapping pulse intervals of width At with the 



state (a, b, c) . These states are distinct, meaning 
a single pulse window At can only be 'in' or 'de- 
scribed by' a single state. Such propositions define 
a state space, which unifies the total time process 
with the total set of configurations. In this scheme 
the total exposure time is decomposed into pulse 
windows of width At, and since every window has 
a state, the temporal composition is equivalent to 
a superposition of independent states. The fraction 
of windows with order k pileup is 



(XAr) k aAt 
k\ 



k k- 



2^Pr((*,fc 

i=Q j=0 



(50) 



where the R.H.S. is the sum over the various com- 
binations with a + b + c = k. In the next section 
it will become clear why each order is sub-divided. 
The terms above are treated as expansion coeffi- 
cients with the spectral components appended. 

The total time process is just the superposition 
of independent states of order k, and thus 



Y (AAT) fc c _ AAT 



fc=0 



k\ 



^Pr(o,6,c) = l (51) 



For example, if there are N events incident dur- 
ing a long exposure time, the average number of 
non-overlapping windows At with configuration 
(tia, tib, Tic) is equal to N * Pr(n A , Ub, Tic)- 



9. Full correction expansion 

The final correction is a combination of the P 
and Q terms and their associated state probabili- 
ties. We first describe the total observed spectrum 
(i.e. the total process), as a superposition of in- 
dependent (non-overlapping) pulse measurements. 
We then devise a correction term based on the fact 
that measurement states overlap, due to the possi- 
bility of recording a tail event. 

9.1. Independent states 

Adopting the assumptions that pileup states are 
independent, the total spectrum can be written as 
a superposition of the peak an d ta il components 

For k th order 



derived given in sections 
pileup, a + b + c = k 



6.3 



and 



7.2 



and the k events can be 
combined into A, B, C to give the set of possible 
pileup states fife = { (a, 6, c) : a + b + c = k}. The 
total number of states per order k is ( fc + 1 )( fc + 2 ) _ 
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The k th order spectrum is an expansion into the 
spectral components with the state probability as 
coefficients. However, since there are k + 1 events 
and a maximum of two can be recorded (one peak, 
one tail) , each p, q term is multiplied by 
the k th order term is written 



1 

k + 1 ' 



Now 



k k—i 



f (k \e) = 



i=Q j=Q 



(52) 



(53) 



where 



p ( o)(e) = S(e) 
Q{a,b,o)( £ ) = 
and the p, q terms have been collected: 

k k—i 

/i%)=EE Pr ^< A ^>*( e ) 

i=0 j=0 
k k—i 

i=0 j=0 



(54) 
(55) 



The total spectrum (PDF), assuming indepen- 
dent pulse states, up to order n, is 

n 

/(e) (56 ) 

= E^l{/i fc) (-) + /f ) ( £ )} (57) 

The approximation order ?i is the value at which 
(A^r)_ e -AAT Decornes negligible. 

9.2. Overlapping pulses (dependent states) 

This expansion is accurate at low rates, when the 
probability of having two adjacent pileup states is 
small. But we must account for the fact that a tail 
measurement causes a new pulse interval to 'begin' 
before the end of the first. A theoretically correct 
construction from the partitioned light curve (i.e., 
the total process) becomes rather complex due to 
this fact. Each state with a tail count, by defini- 
tion, also contains the 'beginning' of the next pulse 




Figure 11: Sample output V(t) showing how exposure time 
is partitioned into pulse states, with measured energies dis- 
tributed as Pa(e) & Qabci^)- Subsequent pulse heights within 
At of a tail measurements are distributed as g a t c (e). A com- 
bination of independent states must be adjusted for the fact 
of overlapping regions. 




Energy [keV] 



Figure 12: Comparison of the spectrum assuming non- 
overalapping pulse states (independent) vs. the spectrum 
corrected for overlap (dependent), at a higher intensity. 
In this case the spectrum is a cutoff power-law, S(E) ~ 
E 15 exp(— E/l MeV). Sigma residuals are calculated rela- 
tive to the monte carlo, assuming Poisson statistics. 



interval. Events occurring within At of the first 
tail count must also be modeled. In terms of the 
state spectral components p and q, the total mea- 
surement from this joined-pulse interval contains a 
peak count, a tail count, and a possible third tail 
count (from the adjoined pulse). 

One way to model this is to consider all possi- 
ble combinations of states. If f2 is the set of pos- 
sible states, then we would have to partition the 
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process in terms of the product states ® f2. At 
very high rates it may even be necessary to consider 
f2®Q®fi, since additional pulse can occur with At 
of the third tail measurement. Due to the obvious 
complications, we have developed a semi-emprical 
approximation which adjust the weighting of peak 
and tail components. 

In the approximation two overlapping intervals 
can be represented as dependent random events M 
and M' , whose probability is generally constructed 
as 

Pr(MUM') =Pr(M)+Pr(M')-Pr(MnM') (58) 

In general the adjoined state, which we may iden- 
tify as M' = (a' , b' , d) can be begin at any random 
time in the interval re, if M records a tail pulse 
M = (a,b,c) with c > 1 and any a,b). Figure 
[ITJ shows several pulse windows with tail measure- 
ments. For a given state (a, 6, c) with c > 1, the 
counts c constitute the peak of a subsequent pileup 
event (c— 1, 6', d). Then we approximate the prob- 
ability Pr(M n M') as the probability of having c 
counts in the tail of a state M and c — 1 counts in 
the peak of a state M' (without regard to the a, b 
or b',d) 



Pr (abc) n ((c- l)b'd) 



(Ar c ) c e^ AT ^ (Ar A ) c - 1 e- ArA 



(c-1)! 



Pr(c, a')5 c ~i. a > 



(59) 



(60) 



The corrected peak contribution is reduced by 
one count per configuration of this tail-peak over- 
lap, so we have a single term subtracted at each 
order from f p 

,(fc) 



D • 

f W (f) - 

J p,Corr\° I 

^{/«( e )_Pr(* + l,fc)p fc (e)} (61) 



k + 



with 



Pr( fc+ l ifc ) = ^C^X^^ (62) 



(k + iy. 



(k) 

The tail correction adds weight to the fq con- 
tribution to the total spectrum, since the adjoined 
pulse has its own tail region. The energy of these 
counts has a similar distribution as the tail of the 



initial pulse, since its baseline is recovered by the 
end of the overlap region (by definition). So we 
model this probability as the probability of having 
a count not in the overlap, over all overlap configu- 
rations; i.e., 1 — Y^k=o P r (^ + 1; k). Thus the total 
'tail' component is 

fq~Corr( £ ) = 

[l - £ Pr(fc' + 1, k')] + / f }(63) 



fc'=0 



fc + 1 



Pr(fc' + l,fc') 



k'=0 



f (fc) 



(64) 



Then the full spectrum is written as 



n 1 



fe=0 



|/W(e)-Pr(fc + l,fc)p fc (e) 



2- £pr(fc' + 1,^)1 



k'=0 



(65) 



or simply 



71 

= E{4So„( £ ) + 4Sorr(^)} ( 66 ) 



fe=0 



This method is evaluated by comparing com- 
puted spectra with monte carlo simulations. Fig- 
ure 12 demonstrates a spectrum computed under 



the assumption of independent pulse states, versus 
the overlap correction. The monte carlo is described 
with more comparisons in section [10.3| Losses due 
to pilcup and subtraction effects can be summarized 
as follows. If the true detection rate is Ao, the pre- 
dicted counting rate (i.e., the rate of pulses heights 
actually being recorded) is 



(67) 



Afl ec = A / f(e)de 

n „ 

= Ao£ / f (k \e)de 

J n " 



fc = ' 



The relationship between Xn ec and Ao is discussed 
further in section [Til 
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10. Numerical evaluation 

In this section we apply the model and compare 
it with monte carlo simulations. The first step is 
numerical evaluation of the probability likelihoods, 
which are used to derive pileup spectral compo- 
nents. 

10.1. Peak time 

The single-pulse peak time is given by f'(tp) = 0. 
Recall that the first-order peak time, defined in gen- 



eral by equation (21), cannot be found in closed 



form for the true pulse shape f(t). Instead we use 
an empirical expression with three constant param- 
eters A s ,A t ,A v to model the dependence of t p on 
two input pulse heights and their separation. 



t 1 p (s 1 ,v (h v 1 )=% + A t - (l- e A^) e - A »^) 2 

(69) 

where t° is the weighted average of the zeroth-order 
peak time of each pulse: 



Vl(t° p + Sl ) 



vq + Vl 



(70) 



and the parameters are manually determined to be 
A„ = |,A S = (|) 2 ,A t = 0.01 fis. Of course these 
parameters, and the functional form itself, are de- 
pendent on the specific pulse shape. The ones given 
here correspond to the best-fit pulse shape param- 
eters from section [3j and 7 is one of the shaping 
parameters from that section. 



10.2. Calculation of Pr^ ak , Pr K A+c , x , B+c 

Because we are using the more accurate pulse 
shape, the conditional probability kernels cannot 
be calculated in closed form. We calculate them 
numerically by first first defining a discrete set of 
channel energies in which probabilities are calcu- 
lated, E = {Ei}. Since pileup events have been 
detected by the instrument it's sufficient to use the 
channel definitions of actual data. 

Defining a discrete time step As, integrals like 



» 



,Pr' 



.<«> 



Pv(e\E,E')= / f sln (s')ds' 



(71) 



are evaluated discretely over the finite As sam- 
ple (s = {0, As,2As, . . . JAs}). At each step of 
the integrand, the function ei(s\E, E') is evaluated. 
Then E\ corresponds to one of the channels in E, 



and an appropriate lookup algorithm returns its in- 
dex (channel) i. For E,E' also discretely sampled 
from E and corresponding to channels j, fc, the con- 
ditional probability is a triply-indexed discrete ob- 
ject Pv(e\E,E') Pr[i,j, k]. 

This probability 'array' is initialized to for all 
i, j, k, and incremented by a value Ap for each step 
in the RHS numerical integration, where I specifies 
the integration step s; to s; + As, and i is the chan- 
nel corresponding to £i(s/| E^, E*,). The interval 
distributions are integrable, so AP/ is exact for the 
I th time step: 



81+1 



AP; 



-(t-s')""W 



= ( T - si) n - (t - s l + 1 ) n 

where r is the normalization interval (A, B, or C) 
and si <E [0,t]. The calculation and storage of the 
Pr(e\E,E') likelihoods introduces some computa- 
tional overhead, since we must have one per order 
in each interval. But since no information about the 
source spectrum or intensity is required, they can 
be calculated a priori and stored in disk memory. 
A separate code computing spectrum corrections 
can read each data block as necessary. In our im- 
plementation we use 128 energy channels and and 
store each function in its own 3D data set in an 
HDF5 file [12]. Using 4-byte floating point data, 
the total uncompressed requirement for all orders 
up to 5 is (3* 128 3 ) x (4 bytes) x 5 = 125MB, which 
is well within the available memory in most envi- 
ronments. 

10.3. Monte carlo comparison 

We have devised a monte carlo to simulate both 
the arrival of random events and the discrete pulse- 
height measurement done by the instrument. We 
simulate a sequence of exponentially distributed 
event times in an arbitrary exposure interval, dur- 
ing which the process intensity is constant. We 
calculate the model up to order 5. For k > 5, a 
trunctation approximation is used by substituting 
fifth-order likelihoods into the iterative calculation. 
At high orders this results in a slight overestimation 
for states having a + b > 5 ^> c. 

An instructive case is that of a narrow line spec- 
trum as might be observed in an ideal laboratory 
experiment. We simulate a gaussian-shaped line 
with fi — 2.2MeV and a = 0.1/^, and compare it 
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with the model prediction (figure 13 ). At low rates 



1000 2000 3000 4000 5000 6000 7000 8000 



the amount of spectral distortion is quite small, but 
as the rate increases the line becomes distorted. Be- 
cause p a (e) and q a bc(s) components are indepen- 
dent of the rate (only their expansion coefficients 
vary), a range of input rates can be tested for a 
given S(E) without much difficulty. 

Plot residuals are calculated in terms of the Pois- 
son variance in each channel of the monte carlo 
output, of = rii. The model is reasonably accu- 
rate for this spectrum, with the average deviation 
< 3<7 in the worst case. For spectral shapes S(E) 
dominated by narrow-band features, such as figure 
13 where the line occupies about 1.6% of the total 
voltage range, the model error tends to be higher 
in the channels away from the line, where the only 
counts are due to high order pileup. To the left 
of the line are primarily tail distorted counts, and 
to the right are mainly counts suffering from peak 
pileup. However for broad-band models, such as the 



one in figure 14 the error is much smaller and more 
uniform. The notable difference is that more chan- 
nels have both peak and tail pileup counts when the 
input spectrum is broad, suggesting that the peak 
and tail modeling assumptions have complementary 
errors, which tend to cancel out when combined. 




2000 3000 
Energy [keV] 



Figure 13: Model vs. simulation comparison for several input 
rates. The input spectrum is a single gaussian-shaped line 
at 2.2 MeV, and is shown by the black dashed line. At high 
rates the line is completely distorted. 
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Figure 14: A spectrum of two lines on a cut-off power-law 
continuum. The model becomes more accurate for contin- 
uum models, suggesting that energy raising and lowering as- 
sociated with peak and tail modeling have canceling errors. 
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Figure 15: A cut-off power law model, whose PDF is ~ 
_E -0 - 5 exp(— _E/7MeV), shown at several rates. The spec- 
trum is plotted on a log-log scale. 



11. Pileup losses 

As already mentioned in the opening sections 
pulse-pileup introduces counting losses which ex- 
ceed the expectation from the conventional non- 
paralyzed deadtime correction. That correction 
is an actually an inversion of the simplest kind, 
one which is algebraically given from the relation 
between the recorded rate and the detected rate, 
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Figure 16: A cut-off power faw modef, whose PDF is ~ 
E 1 - 5 exp(— E/lMeV), shown at severai rates on a log-log 
scale. 
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Figure 17: A spectrum of four fictitious emission lines, with 
no continuum. Despite large error in the 38 kcps spectrum, 
the predicted number agrees with simulation. 



Afl ec (Ao), in the presence of a non-paralyzable dead- 
time r. By simple arguments this function is Xn ec — 
Aq/(1 + Aqt), and is algebraically invertible to give 



Ao — A_R ec /(l — AflecT") 



(72) 



At high rates it is insufficient both because of 
pileup, and because the derivation assumes Xr &c < 



Figure 18: A corrected spectrum of emission lines on an 
exponential continuum, compared with the monte carlo sim- 
ulation at several different input rates. The first three order 
terms /W (e), /( 2 ) (e), /( 3 ) (s) with the rate-dependent nor- 
malization arc also shown. 



i/r. n, m 

Pileup changes the picture for two reasons. The 
first is that it randomly extends the instrument 
deadtime. In GBM this can be considered as a kind 
of paralyzable deadtime since the additional detec- 
tions delay peak measurement. The second is that 
the tail effect causes energy-dependent baseline sub- 
traction losses which are more significant for 'flat' 
or broadband spectra having a mixture of large and 
small pulses. These effects must be accounted for 
to fully correct the relation Aij ec (An). 

The first effect (paralyzable deadtime) by itself 
can be predicted, to first order, using Poisson statis- 
tics and employing a root-finding algorithm. The 
procedure is described in [5] , and yields a numerical 
solution. We do not explicitly use this treatment in 
our model. Rather, because the model's predicted 
losses agree well with simulation, we conclude that 



the method of overlapping windows, equation (65), 



is an effective proxy for paralyzable deadtime. This 
seems reasonable since the subtracted overlap terms 



in equation ( 6f I reduce the proportion of peak mea- 



surements, which is roughly the same effect as par- 
alyzeble deadtime. An alternative interpretation is 
that severe paralyzable deadtime is improbable in 
GBM even at the high rates tested, due to the pulse 
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shape and peak finding algorithm. At rates beyond 
10 6 cps, it is likely that the instrument becomes 
non-linear and our assumptions would not hold. 

The second effect, tail subtraction loss, is some- 
what more serious than the first in our case, be- 
cause the pulse tail is longer than the peak and the 
negative amplitude is large. Because of spectral 
dependence, tail losses can only be predicted by 
assuming a zeroth order pulse-height distribution 
S(E). These losses can be investigated by plotting 
marginal distributions from the tail likelihoods cal- 
culated in section [7] Figures 21 and 22 show the 
first-order case of the 'A+C and ! B+C pileup sce- 
narios. When tail spectral components are calcu- 



lated using probability integrals like equation (47), 



they turn out to have total probability less than 
one (J2 E Q{s) < 1). Evidently this due to regions 
of where losses are > 0%. By contrast the peak 
components are all normalized to one. 

The model we present gives an accurate predic- 
tion for the additional pulse-pileup losses. Figure 



19 shows that the model result from equation (65) 
is consistent with monte carlo simulation. Figure 

20 shows Aij ec (Ao) for the full range of tail losses 
due to variations in spectral hardness. Most real 
spectra will be in the upper half of the shaded re- 
gion. 

At present an analytical inversion giving Ao(Afl ec ) 
with deadtime and pileup has not been found, 
though in principle one exists until the the turnover 
in figure |20| Future efforts using a technique such 



as series reversion of equation ( 65 ) may be fruit 



ful. However we note that spectral fitting with the 
pileup correction is itself an inversion process and 
results in a possible solution for Ao given Xn ec from 
real data. 




Figure 19: Fraction recorded vs. detection rate. 
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12. Conclussions 

We have derived a model which accurately pre- 
dicts the recorded number and spectrum of a con- 
stant intensity Poisson process with pulse-pileup, 
when compared to monte carlo simulations. We 
have used the peak modeling technique and iter- 
ation method of [I], and extended the treatment 
to a three-region bipolar pulse. The novelty in this 
method is that it provides a way to model input en- 
ergy and timing statistics of tail pileup events. The 
total spectrum is written as a state-space expansion 
of overlapping pulses. The technique generally ap- 
plies to bipolar shaping instruments, and has been 
demonstrated using the true pulse shape for GBM. 
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Figure 20: Relation between the detection rate Ao and the 
observed rate \R ec , in millions of counts per second. As the 
detection rate increases, a simple deadtime approximation 
becomes insufficient to explain the recorded rate. Moreover, 
energy dependent losses occur due to the tail effect, which 
are worse for broadband, energetic ("flat") spectra because 
they generate small and large pulses. Pulse pileup greatly 
increases uncertainty about the true rate given the observed 
rate. As the input rate exceeds 2* 10 cps (Aot 0.5), con- 
straining Ao given \n ec becomes difficult (only a lower limit 
can be determined). Note, since GBM detectors have their 
own counting electronics, Aq is the rate in a single detector. 
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Figure 21: Baseline subtraction losses due to the (first or- 
der) A+C effect. Eq is energy of the primary count. E\ is 
the input energy of the second count, which occurs in the 
tail of Eq (i.e., region C). This is the same configuration as 
depicted in figures [5Jc) and [9] As the peak amplitude Eq 
increases, its tail becomes more negative. If E\ is too small, 
the summed peak in C is below registration threshold and E\ 
is lost. 100% loss occurs when Eq 2> E\. Each curve is cal- 
culated directly from the tail likelihood by integrating over 
the recorded pulse height (sum over channels in the discrete 
case): Pr loas (E , E^ = (1 - £ Pr A +C^\E , £i)) 
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Figure 22: Losses in the B+C effect (figure |10) . Each 
curve is calculated by integrating the measurement prob- 
ability in the event of Eq,Ei: Pri oss {Eq, Ei) = (1 — 
E^ s+ c(£|^0,Si)). 
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